library(facetscales)


#effect of failure analyses------------------------------------------------------------------------------------------------------------------------------------------
estimation_sample_javanese = data_1 %>% filter(abs(forcing_skd) < 5, java_indicator == 1)
estimation_sample_non_javanese = data_1 %>% filter(abs(forcing_skd) < 5, java_indicator == 0)
estimation_sample = data_1 %>% filter(abs(forcing_skd) < 5)

java1_java = lm_robust(q17_survey_java1 ~ fail_skd, data = estimation_sample_javanese) %>% tidy() %>% mutate(outcome = "Govt should focus attention on Java (1-4)")
java2_java = lm_robust(q17_survey_java2~ fail_skd, data = estimation_sample_javanese) %>% tidy() %>% mutate(outcome = "Govt has given most resources to Java (1-4)")

java1_nonjava = lm_robust(q17_survey_java1~ fail_skd, data = estimation_sample_non_javanese) %>% tidy() %>% mutate(outcome = "Govt should focus attention on Java (non, 1-4)")
java2_nonjava = lm_robust(q17_survey_java2 ~ fail_skd, data = estimation_sample_non_javanese) %>% tidy() %>% mutate(outcome = "Govt has given most resources to Java (non, 1-4)")

region1 = lm_robust(q18_survey_daerah1 ~ fail_skd, data = estimation_sample) %>% tidy() %>% mutate(outcome = "Local govt shoud focus attention on locals over immigrants (1-4)")
region2 = lm_robust(q18_survey_daerah2 ~ fail_skd, data = estimation_sample) %>% tidy() %>% mutate(outcome = "Too many outsiders work in government (1-4)")
region3 = lm_robust(q18_survey_daerah3 ~ fail_skd, data = estimation_sample) %>% tidy() %>% mutate(outcome = "Local govt focuses too much on city-dwellers (1-4)")

relg1 = lm_robust(q19_survey_relg1  ~ fail_skd, data = estimation_sample) %>% tidy() %>% mutate(outcome = "Upset if different religion built place of worship (1-4)")
relg2 = lm_robust(q19_survey_relg2  ~ fail_skd, data = estimation_sample) %>% tidy() %>% mutate(outcome = "Upset if different religion became mayor (1-4)")
relg3 = lm_robust(q19_survey_relg3  ~ fail_skd, data = estimation_sample) %>% tidy() %>% mutate(outcome = "Upset if different religion became national minister (1-4)")

nation1 = lm_robust(q20_survey_pancasila_h3 ~ fail_skd, data = estimation_sample) %>% tidy() %>% mutate(outcome = "Is Pancasila still relevant? (1-4)")
nation2 = lm_robust(q21_survey_nation_ethnic_h3 ~ fail_skd, data = estimation_sample) %>% tidy()%>% mutate(outcome = "Identify more as ethnic group (1) or Indonesian (3)")

corrupt1 = lm_robust(q16_survey_transparent_h2  ~ fail_skd, data = estimation_sample) %>% tidy()%>% mutate(outcome = "How transparent is recruitment? (1-4, reversed)")
corrupt2 = lm_robust(q13_survey_success_reason_1_merit_h2 ~ fail_skd, data = estimation_sample) %>% tidy()%>% mutate(outcome = "How important is merit (1-4, reversed)")
corrupt3 = lm_robust(q13_survey_success_reason_3_sara  ~ fail_skd, data = estimation_sample) %>% tidy()%>% mutate(outcome = "How important is ethnicity (1-4)")
corrupt4 = lm_robust(q13_survey_success_reason_2_connection_h2  ~ fail_skd, data = estimation_sample) %>% tidy()%>% mutate(outcome = "How important are connections (1-4)")
corrupt5 = lm_robust(q14_survey_test_vs_connection ~ fail_skd, data = estimation_sample) %>% tidy() %>% mutate(outcome = "Which more important? Test (0) or connections (1)")



#service effect analyses------------------------------------------------------------------------------------------------------------------------------------------

estimation_sample_javanese = data_2 %>% filter(abs(forcing_final) < 1, java_indicator == 1)
estimation_sample_non_javanese = data_2 %>% filter(abs(forcing_final) < 1, java_indicator == 0)
estimation_sample = data_2 %>% filter(abs(forcing_final) < 1)

java1_java_s = lm_robust(q17_survey_java1 ~ pass_final, data = estimation_sample_javanese) %>% tidy() %>% mutate(outcome = "Govt should focus attention on Java (1-4)")
java2_java_s = lm_robust(q17_survey_java2~ pass_final, data = estimation_sample_javanese) %>% tidy() %>% mutate(outcome = "Govt has given most resources to Java (1-4)")

java1_nonjava_s = lm_robust(q17_survey_java1~ pass_final, data = estimation_sample_non_javanese) %>% tidy() %>% mutate(outcome = "Govt should focus attention on Java (non, 1-4)")
java2_nonjava_s = lm_robust(q17_survey_java2 ~ pass_final, data = estimation_sample_non_javanese) %>% tidy() %>% mutate(outcome = "Govt has given most resources to Java (non, 1-4)")

region1_s = lm_robust(q18_survey_daerah1 ~ pass_final, data = estimation_sample) %>% tidy() %>% mutate(outcome = "Local govt shoud focus attention on locals over immigrants (1-4)")
region2_s = lm_robust(q18_survey_daerah2 ~ pass_final, data = estimation_sample) %>% tidy() %>% mutate(outcome = "Too many outsiders work in government (1-4)")
region3_s = lm_robust(q18_survey_daerah3 ~ pass_final, data = estimation_sample) %>% tidy() %>% mutate(outcome = "Local govt focuses too much on city-dwellers (1-4)")

relg1_s = lm_robust(q19_survey_relg1  ~ pass_final, data = estimation_sample) %>% tidy() %>% mutate(outcome = "Upset if different religion built place of worship (1-4)")
relg2_s = lm_robust(q19_survey_relg2  ~ pass_final, data = estimation_sample) %>% tidy() %>% mutate(outcome = "Upset if different religion became mayor (1-4)")
relg3_s = lm_robust(q19_survey_relg3  ~ pass_final, data = estimation_sample) %>% tidy() %>% mutate(outcome = "Upset if different religion became national minister (1-4)")

nation1_s = lm_robust(q20_survey_pancasila_h3 ~ pass_final, data = estimation_sample) %>% tidy() %>% mutate(outcome = "Is Pancasila still relevant? (1-4)")
nation2_s = lm_robust(q21_survey_nation_ethnic_h3 ~ pass_final, data = estimation_sample) %>% tidy()%>% mutate(outcome = "Identify more as ethnic group (1) or Indonesian (3)")

corrupt1_s = lm_robust(q16_survey_transparent_h2  ~ pass_final, data = estimation_sample) %>% tidy()%>% mutate(outcome = "How transparent is recruitment? (1-4, reversed)")
corrupt2_s = lm_robust(q13_survey_success_reason_1_merit_h2 ~ pass_final, data = estimation_sample) %>% tidy()%>% mutate(outcome = "How important is merit (1-4, reversed)")
corrupt3_s = lm_robust(q13_survey_success_reason_3_sara  ~ pass_final, data = estimation_sample) %>% tidy()%>% mutate(outcome = "How important is ethnicity (1-4)")
corrupt4_s = lm_robust(q13_survey_success_reason_2_connection_h2  ~ pass_final, data = estimation_sample) %>% tidy()%>% mutate(outcome = "How important are connections (1-4)")
corrupt5_s = lm_robust(q14_survey_test_vs_connection ~ pass_final, data = estimation_sample) %>% tidy() %>% mutate(outcome = "Which more important? Test (0) or connections (1)")



bold_labels <-  c("Javan Preferentialism (Javans):                                          ", "Javan Preferentialism (non-Javans):                                      ",
                  "Regional Preferentialism:                                                   ", "Religious Intolerance:                                                      ",
                  "National Identification:                                                    ", "Perceptions of Corruption:                                                  ")

plot_data_main <-
  bind_rows(java1_java, java2_java, java1_nonjava, java2_nonjava,region1, region2, region3, 
            relg1, relg2, relg3, nation1, nation2, corrupt1, corrupt2, corrupt3, corrupt4, corrupt5) %>%
  filter(term == "fail_skd") %>%
  left_join(.,
            bind_rows(java1_java, java2_java, java1_nonjava, java2_nonjava,region1, region2, region3, 
                      relg1, relg2, relg3, nation1, nation2, corrupt1, corrupt2, corrupt3, corrupt4, corrupt5) %>%
              filter(term == "(Intercept)") %>%
              dplyr::select(outcome, cntrl_mean = estimate)) %>%
  dplyr::select(estimate, std.error, outcome, p.value, cntrl_mean) %>%
  mutate(class = "Effect of Failure")


plot_data_s <-
  bind_rows(java1_java_s, java2_java_s, java1_nonjava_s, java2_nonjava_s,region1_s, region2_s, region3_s, 
            relg1_s, relg2_s, relg3_s, nation1_s, nation2_s, corrupt1_s, corrupt2_s, corrupt3_s, corrupt4_s, corrupt5_s) %>%
  filter(term == "pass_final") %>%
  left_join(.,
            bind_rows(java1_java_s, java2_java_s, java1_nonjava_s, java2_nonjava_s,region1_s, region2_s, region3_s, 
                      relg1_s, relg2_s, relg3_s, nation1_s, nation2_s, corrupt1_s, corrupt2_s, corrupt3_s, corrupt4_s, corrupt5_s) %>%
              filter(term == "(Intercept)") %>%
              dplyr::select(outcome, cntrl_mean = estimate)) %>%
  dplyr::select(estimate, std.error, outcome, p.value, cntrl_mean) %>%
  mutate(class = "Effect of Success")


plot_data <-
  bind_rows(plot_data_main, plot_data_s) %>%
  add_row(., estimate = NA_real_, std.error = 0, p.value = 1, class = "Effect of Success", outcome = "Javan Preferentialism (Javans):                                          ") %>%
  add_row(., estimate = NA_real_, std.error = 0, p.value = 1, class = "Effect of Success", outcome = "Javan Preferentialism (non-Javans):                                      ") %>%
  add_row(., estimate = NA_real_, std.error = 0, p.value = 1, class = "Effect of Success", outcome = "Regional Preferentialism:                                                   ") %>%
  add_row(., estimate = NA_real_, std.error = 0, p.value = 1, class = "Effect of Success", outcome = "Religious Intolerance:                                                      ") %>%
  add_row(., estimate = NA_real_, std.error = 0, p.value = 1, class = "Effect of Success", outcome = "National Identification:                                                    ") %>%
  add_row(., estimate = NA_real_, std.error = 0, p.value = 1, class = "Effect of Success", outcome = "Perceptions of Corruption:                                                  ") %>%
  mutate(outcome = factor(outcome, levels = 
                            c("Javan Preferentialism (Javans):                                          ", "Govt should focus attention on Java (1-4)", "Govt has given most resources to Java (1-4)",
                              "Javan Preferentialism (non-Javans):                                      ", "Govt should focus attention on Java (non, 1-4)", "Govt has given most resources to Java (non, 1-4)",
                              "Regional Preferentialism:                                                   ", "Local govt shoud focus attention on locals over immigrants (1-4)", "Too many outsiders work in government (1-4)", "Local govt focuses too much on city-dwellers (1-4)",
                              "Religious Intolerance:                                                      ", "Upset if different religion built place of worship (1-4)", "Upset if different religion became mayor (1-4)", "Upset if different religion became national minister (1-4)",
                              "National Identification:                                                    ", "Is Pancasila still relevant? (1-4)", "Identify more as ethnic group (1) or Indonesian (3)",
                              "Perceptions of Corruption:                                                  ", "How transparent is recruitment? (1-4, reversed)", "How important is merit (1-4, reversed)", "How important is ethnicity (1-4)", "How important are connections (1-4)", "Which more important? Test (0) or connections (1)"
                            ))) %>%
  mutate(outcome = fct_rev(outcome))

bold.labels <- ifelse(levels(plot_data$outcome) %in% bold_labels, yes = "bold", no = "plain")

scales_x <- list(
  `Effect of Failure` = scale_x_continuous(limits = c(-0.15, 0.25), breaks = seq(-0.15, 0.25, .1)),
  `Effect of Success` = scale_x_continuous(limits = c(-0.6, 0.3), breaks = seq(-0.6, 0.3, 0.15))
)

individual_items <-
  plot_data %>%
  mutate(label = case_when(p.value < 0.1 & estimate > 0 ~ paste0(round(estimate, 3), " (+", round((estimate/cntrl_mean)*100, 1),"%", ")"),
                           p.value < 0.1 & estimate < 0 ~ paste0(round(estimate, 3), " (", round((estimate/cntrl_mean)*100, 1),"%", ")"),
                           TRUE ~ " ")) %>%
  ggplot(aes(y=outcome, x = estimate)) +
  geom_point() +
  geom_text(aes(x = case_when(estimate > 0  & class == "Effect of Failure" ~ (estimate+1.67*std.error) + 0.05,
                              estimate > 0  & class == "Effect of Success" ~ (estimate+1.67*std.error) + 0.1,
                              estimate < 0  & class == "Effect of Failure" ~ (estimate-1.67*std.error) - 0.05,
                              estimate < 0  & class == "Effect of Success" ~ (estimate-1.67*std.error) - 0.1,
  ),
  label = label), 
  size = 3) +
  geom_errorbarh(aes(xmin = estimate - 1.68*std.error, xmax = estimate + 1.68*std.error), size = 0.8, height = 0) +
  geom_errorbarh(aes(xmin = estimate - 1.96*std.error, xmax = estimate + 1.96*std.error), size = 0.4, height = 0) +
  geom_vline(xintercept = 0) +
  theme_bw() +
  facet_grid_sc(cols = vars(class), scales = list(x = scales_x)) +
  theme(panel.grid.minor = element_blank(), 
        panel.grid.major.y = element_blank(),
        axis.line.y.left = element_blank(),
        axis.text.y = element_text(face = bold.labels),
        axis.title.x = element_blank(),
        axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        legend.title = element_blank(),
        axis.ticks.y = element_blank(),
        strip.background = element_blank(),
        legend.position = "none") +
  theme(text=element_text(size=11)) +
  ylab("")

ggsave(filename = "./_4_outputs/figures/figure_3.pdf", plot = individual_items, width = 13, height = 6.5)
